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ABSTRACT 

We write down and solve equations describing steady state, optically thin, 
advection-dominated accretion onto a Kerr black hole. The mean flow, described 
by the relativistic fluid equations, is axisymmetric and vertically averaged. The 
effect of turbulence in the flow is represented by a viscous shear stress. Our 
treatment differs in several important ways from earlier work: we use a causal 
prescription for the shear stress, we do not assume the relativistic enthalpy is 
unity (this is important for rapidly rotating holes), and we use a relativistic 
Q\ . equation of state. We present several representative solutions and use them 

to evaluate the importance of relativistic effects, to check our approximations, 
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and to evaluate the robustness of the input physics. Detailed properties of the 
solutions are described in an accompanying paper. 



1. Introduction 

The structure of accretion flows close to the horizon of a black hole is of considerable 
astrophysical interest because: (1) most of the gravitational binding energy is released close 
to the black hole; (2) there are strong-field gravitational effects close to the horizon that are 
unique to black holes; (3) temperatures are highest close to the horizon, so exotic physical 
processes may occur there. 



Early studies of disk accretion onto black holes assumed a thin disk. |Shakura~£; 



Sunyaev ( 1973 )| presented the basic equations used to describe thin disks, but did not 



include relativity. Accretion in thin, relativistic disks was first described by Pardeen, Press 



fe Teukolsky (1972)| and [INovikov fc Thome (19731 - Note that recently |Rittert fc Herold 



(1995) corrected an algebraic error in NT's version of the vertical structure equation. 
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Paczyhski & Bisnovatyi-Kogan (1981), [Muchotrzcb fc Paczynski (1982) , and 



Abramowicz et al. (1988)] introduced a set of disk equations, dubbed the "slim disk" 



equations by |Abramowicz et al. (1988)1 , which include a number of terms neglected in the 



thin disk formulation. These equations include radial pressure gradient, acceleration, and 
energy transport, and provide the basis for the equations used in this paper and most of the 
other recent work cited here on disk structure around black holes. These studies all used 
the "Paczynski potential", a pseudo-relativistic potential designed to reproduce the main 
features of orbits around a nonrotating black hole ( Paczynski fc Wiita (1980~)| ). 

One of the terms included in the slim disk equations is the inward advection of 
entropy by the accreting gas. [Narayan fc Popham (199UJ] showed that this term, which is 



insignificantly small in the usual thin disk, becomes large in the optically thin boundary 
layer region of low-accretion rate cataclysmic variables. Narayan and Yi recognized the 
existence of a class of flows in which advection dominates the energy balance, and found 
a self-similar solution in the advection-dominated limit ( [Narayan fc Yi (1994)| , [Narayan 



|fc Yi (1995a)| ). Advection-dominated solutions were recognized as a new class of thermal 
equilibria QAbramowicz et al. (1995) , Uhen et al. {1995) , Narayan fc Yi (1995b) ). They 



resemble some earlier models by |Rees et al. (1982)] and |Bcgclman fc Meier (1982~| , where 



the hot accreting gas radiates very inefficiently. 

In this paper we concentrate on a particular class of flows: optically thin, advection- 
dominated accretion flows, or ADAFs. These flows cool slowly, so most of the heat deposited 
in them by the dissipation of turbulence (which also carries off angular momentum) is 
advected into the black hole rather than radiated away. The accretion is therefore inefficient 
compared to a thin disk, which can radiate away a substantial fraction of the accreted rest 
mass energy. Over the past few years, the theory of these flows has been developed and 
applied to a number of observed black hole systems (see |JNarayan (1997)| for a review). Most 



recently, |Narayan, Kato, fc Honma (1997)| (NKH) and |Chen, Abramowicz, fc Lasota (1997) 



presented global, advection-dominated disk solutions in the Paczynski potential. Our work 
differs from this in that it is fully relativistic, and includes effects due to strong rotation of 
the black hole. 

Recently, a number of authors have included the effects of relativity in their models of 
disks around black holes. [Lasota (1994)| first wrote down slim disk equations which included 



relativistic effects. This was superseded by the work of |Abramowicz et al. (1997)| (ACGL), 



who corrected some errors in Lasota's work and provided sample numerical solutions. Our 
work differs from this one in many details, but mainly in our use of a turbulent stress 
prescription that is explicitly related to the Navier-Stokes stress and that is causal, and in 
our inclusion of the relativistic enthalpy. |Abramowicz, Lanza, fc Percival (1997)| (ALP) 
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revisited the question of vertical structure and derived a new expression for the scale height 
beginning with the Euler equations. We adopt ALP's version of the vertical structure 
equations, with a minor correction. Independently, [Peitz fc Appl (1997)| (PA) gave a very 
nice and rather complete treatment of the problem, which included a causal viscosity and 
the relativistic enthalpy; however, they used a polytropic equation of state rather than 
solving the energy equation as we do. |Chakrabarti (1996)1 has also considered rotating 
accretion flows in the Kerr metric, but only in the weak viscosity limit. Also, unlike all of 
the above, we use a relativistic equation of state. This correction increases in importance 
as the hole rotates more and more rapidly. Finally, we correct some minor errors in some of 
the above treatments. 

The spirit of this work is to do the best possible job on the physics, including the 
relativity and turbulent shear stress, within certain limitations. One limitation is physical: 
there is no well justified formalism for treating the effects of turbulence in these flows. We 
treat the mean flow using the full equations of relativistic hydrodynamics, and include 
the effects of turbulence via a turbulent shear stress, making the minimal modifications 
of the inviscid equations necessary to allow angular momentum transport and therefore 
accretion. Another limitation is complexity: we are forced to make some approximations 
simply because without them the equations would be too difficult to solve. Despite these 
limitations, this work reveals some generic properties of relativistic ADAFs, and, when 
coupled with a proper scheme for relativistic photon transport, will allow us to model their 
observational appearance. 

The plan of this paper is as follows. In §2 we explain our units, notation, and record 
the relevant metrics and frames. In §3 we write down the basic dynamical equations. In 
§4 we consider the heart of the physics in this problem, the turbulent shear stress. In §5 
we evaluate the critical point conditions. In §6 we give a sample solution (a full survey of 
solutions is described in |Popham fc Gammie (1997)| ), and evaluate our approximations. §7 
contains our conclusions. 



2. Preliminaries 

2.1. Units and Metric 

We adopt G = M = c= las our basic scalings, where M is the black hole mass. 
This implies a unit of mass, length, and time of M, GM/c 2 , and GM/c 3 , respectively. The 
angular momentum J of the black hole is described by a = Jc/GM 2 , where — 1 < a < 1. 
It is convenient to define, like NT, the following relativistic correction factors that become 
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unity in the nonrelativistic limit: 

.4=1 + a 2 /r 2 + 2a 2 /r 3 , (1) 

C = 1 - 3/r + 2a/r 3 / 2 , (2) 

D = 1 - 2/r + a 2 /r 2 , (3) 

and 

// = 1 + a 2 cos 2 #/r 2 . (4) 
In Boyer-Lindquist coordinates, the Kerr metric is 

ds 2 = -(l-^)dt 2 -^f^dtd<f) + %dr 2 + r 2 iid6 2 + (5) 

r 2 sin 2 #(l + J + 2aW*)d0 2 . (6) 



The non-zero contravariant components are 

= -1-2(1 + a 2 /r 2 )/r£>/i, g^ = -2a/r 3 V^ g rr = V/fi, 

g °6 = l/ r 2 /i, ^ = (1 - 2/r / u)/r 2 sin 2 #£>. 

Often we shall only require the metric in the equatorial plane, where 



(7) 



ds 2 = -^-dt 2 + r 2 A(d<j) - uodtf + ]-dr 2 , (8) 



and 

measures the rate of frame dragging by the hole. In expanded form, 

ds 2 = -(1 - -)dt 2 - —dtd(f) + -dr 2 + r 2 Ad<p 2 . (10) 
r r V 

The nonzero contravariant components of the metric in the equatorial plane are 

g u = -A/V, = -2a/r 3 V, g rr = V, = (1 - 2/r)/r 2 V. (11) 

In the equatorial plane, the horizon lies at the outer root of T> — 0, i.e. r = 1 + y/l — a 2 . 
The boundary of the "ergosphere" , where the world lines of observers at constant 
Boyer-Lindquist coordinates (r, 6, <p) become spacelike (and therefore unphysical) lies at 
r = 2 in the equatorial plane. 
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2.2. Basic Flow Symmetry, Vertical Averaging 

Our goal is to accurately describe the flow of fluid in the highly relativistic regime, close 
to the horizon. To do so we must make some simplifying assumptions. First, we assume 
the angular momentum of the accreting fluid is aligned with the angular momentum of the 
black hole. f\ Second, we assume that the flow, described by its four-velocity (V, u r , u e , uP) 
is in the mean axisymmetric and symmetric about the equatorial plane. Then uq vanishes 
at the equatorial plane. 

We also vertically average the flow, a standard approximation in accretion flows with 
angular momentum. This is best explained by example. Consider the particle number 
conservation equation: 

(pu^ = 0. (12) 

Here p is the "rest-mass density" . Then 

(pit");, = 4=(V^").A« = l(^V),r + — !-r(//sin0 A* = 0, (13) 

where g = |Det(^)| = r 4 sin6 |2 /x 2 . 

We now perform a "vertical" average by integrating over the volume between r and 
r + Sr. Then the second term in equation (|T3| ) vanishes. Now define He, the characteristic 
angular scale of the flow about the equator, and assume this to be the same for all flow 
variables /. The vertical averaging approximation consists in taking 

J m<f>y/gS ^ AnHefie = vr/2), (14) 
(4nH e r 2 pu r ), r = 0. (15) 



whence 



Integrating once in radius, 

AnH e r 2 pu r = —M, (16) 
where the constant M is the "rest-mass accretion rate". 



2 We make this assumption because the full non-aligned problem is intractable. Alignment of the flow due 
to gravitomagnetic precession and viscosity (the Bardeen-Petterson effect, see Bardeen fc Petterson (1975)| , 



Kumar fc Pringlc (1985)| ) is not likely to be important for the hot flows considered here. 
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2.3. Frames 

We use four separate frames in our calculations. The first and most important is the 
Boyer-Lindquist coordinate frame (BLF), where most direct calculations are done. 

Second is the locally nonrotating frame (LNRF), an orthonormal tetrad basis carried 
by observers who live at constant 9 and r, but at <p = ut + const., so they are dragged in 
azimuth by the black hole. Their world lines are always timelike, even in the ergosphere. 
Explicit transformations between the LNRF and BLF are given by |Bardeen, Press, 
[leukolsky (19721 

A third frame is obtained by an azimuthal Lorentz boost from the LNRF into another 
tetrad basis that corotates with the fluid (the corotating frame, or CRF). The CRF has 
velocity (3$ with respect to the LNRF; we also define 7^ = (1 — /3|) -1 / 2 . 

Finally we have the local rest frame (LRF) of the fluid, obtained by a radial Lorentz 
boost from the CRF. The LRF has radial velocity V < with respect to the CRF, so we 
define (3 r = V and j r = (1 - /3 r 2 )" 1/2 . 

Two velocity variables are needed to describe the flow in the equatorial plane. We 
choose V (as do ACGL), the radial velocity of the fluid measured in the CRF, and I = u<j> 
(as do PA), the angular momentum. 

All the contravariant and covariant components of the velocity field can now be 
expressed in terms of I and V. We find, using 7 = 7^7,., 



(u\u r ,u\v?) 
In terms of / and V, 

2 _ 1 I 2 
7 ~ 1-V 2 + r 2 A' 

^ r\fAr) 

Recall that £ = —u t is the "energy at infinity", and that 
geodesies. 

Finally, it is useful to have explicit expressions for the LRF basis vectors. The 



V , V 
V A ,/l?(l - V 2 ) 



0, I 



71 



V\ 



I V 



r 2 A 



+ 0JJ\ 



;i7) 

18) 



(19) 
(20) 

/ and £ are conserved along 
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contravariant basis vectors are 



e^ = (Ji, PrlrVV, 0, ^2= + ^=) 

(t) wy p > r^y/AD rVA 



A r^z 2aP r 'j P r 



e M (r) = Wrl\ 7^, IrVV, 0, 



V 



+ 



e>) = (0, 0, -, 0) 



The covariant basis vectors are 



, (*) = {l K+M£L o -B^y/A) 



-4 r 2 ^4 ' VV 
eW = (0, 0, r, 0) 



, 0, PtPr-yrVA) 



(21) 

(22) 
(23) 
(24) 

(25) 

(26) 
(27) 
(28) 



Recall that the basis vectors allow one to transform back and forth from the LRF via 



3. Basic Dynamical Equations 

The basic equations of relativistic, viscous hydrodynamics are the energy-momentum 
conservation equations 

T^, v = 0, (29) 
and the particle number conservation, or continuity, equation: 

(pu 11 );, = 0. (30) 

We adopt the convention that p is the "rest mass density" and u is the internal energy per 
unit proper volume, so the total density of mass-energy is p + u. We also define the pressure 
p and the relativistic enthalpy rj = (p + u + p) / p. Then the stress-energy tensor is 



(31) 
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where is the "viscous" stress-energy tensor. It includes terms due to small-scale velocity 
and magnetic field fluctuations. We neglect contributions to the stress-energy tensor from 
the energy flux and from large-scale electromagnetic fields such as those found in MHD 
winds. 

In relativistic Navier-Stokes flow with zero bulk viscosity 

1T = -2\a^, (32) 
where A is the "coefficient of dynamic viscosity" , f\ and 

Va/3 = - (uawh'p + u^h^j - -Qh a/3 , (33) 

is the shear tensor. Here 

6 = u%, (34) 

and 

is the projection tensor. Our turbulent shear stress tensor does not have this form, that is, 
it is not a relativistic Navier-Stokes viscous stress tensor. Instead we write down a turbulent 
stress tensor, described in §4, that is both the minimal departure from the inviscid equations 
that allows angular momentum transport and the simplest form that preserves causality. 



3.1. Particle Number Conservation 

As already discussed, the particle number conservation equation can be reduced to the 

form 

- 47rr 2 /m r H e = M. (36) 
Expressing u r in terms of the preferred dependent variable V, 

Anr 2 pH e V (^— -j = —M. (37) 

A closely related equation may be derived with the aid of the Killing vector 
# = (1,0,0,0). Then 

= = ~ 2 {r 2 T\), r + — L^sm0T* t V (38) 
r z fi sin 



Israel (1972) 



for a discussion of relativistic microscopic viscosity. 
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Vertically averaging, the second term vanishes. Then integrating once in r and expanding 
we find 

4TrH e r 2 [-(p + u + p)£u r + t t r ] = E. (39) 

This equation expresses the constancy of mass-energy flux E with radius; E is the actual 
rate of change of the black hole mass. If the flow is cold and slow at large radius, E ~ M. 



Subtracting equation (^g) from equation (|39|) then dividing by equation (p6|) gives the 
following relation: 

_ ^ + £ -l + (^l)s=^L. (40) 
pu r V p J M 

This reduces to the Bernoulli equation in the inviscid limit. 



3.2. Energy Equation 

The energy equation is the component of equation (p9| ) parallel to the fluid four- velocity: 

n^ v , v = 0. (41) 
Using standard manipulations (e.g. [Ellis (197T)|) , this can be reduced to 



„-^_ u r(«±rf^ = $ _ A (42) 

dr p dr 

where, provided that t^ v is trace-free, as we shall assume, 

$ = -tTa^ (43) 

is the dissipation function and A is the cooling function. 

Two variables are required to represent the thermodynamic state of the fluid. We 
choose p, the rest-mass density, and the scaled temperature T. In dimensional form, T is 
related to the ordinary kinetic temperature Tk by T = fcT^/mc 2 , where k is Boltzmann's 
constant and fh is the mean molecular weight. 

The pressure is given by the ideal gas equation of state p = pT. We fold the magnetic 
contribution to the pressure into p and any magnetic contribution to the energy per unit 
proper volume into the internal energy u. We must now specify u(p,T). We require an 
equation of state that is polytropic in the nonrelativistic limit with adiabatic index j . 
Since the magnetic energy is included in the internal energy, 7 Q may differ from 5/3. Thus 
we set 

T _ /4/(to-1) + 15T\ 
u = P Tg(T) = pT I 1 . (44) 
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This formula fits the exact relativistic equation of state for an ideal relativistic Boltzmann 
gas (see, e.g., [Service (1986)| ), for which 7 = 5/3, to everywhere better than 2%. [] The 
sound speed is then 

2 Tp TT 

c = — = — , 45 
rjp 7] 

where 



_ dlnp 



1+ , T A JJfTy (46) 

g + T(dg/dT) 



dlnp 

When 7 = 5/3, our T agrees with Service's to within 0.5%, and our sound speed agrees 
with hers to within 1.3%. 

The energy equation can now be reduced to 

T / V \ 1/2 fdudT pdp\ ^ 



. 1 — V 2 J \ dT dr p dr 
The dissipation function can be evaluated once we specify the viscous stress. 

Like NKH, we do not include cooling explicitly, but set $ — A = /$, where / < 1 is 
a constant factor that we vary to estimate the effect of cooling. Our energy equation is 
essentially the same as ACGL, except that we include the relativistic equation of state, and 
our dissipation function is different. PA do not solve an energy equation; instead they use a 
polytropic equation of state. 



3.3. Radial Momentum 

The radial momentum equation is the radial component of the projection of equation 
(fZ§|) into the space normal to the four-velocity: 

M^"");, = 0. (48) 

Using standard manipulations (we must evaluate some connection coefficients), and 
neglecting the viscous acceleration, the radial momentum equation becomes 

V ' 1V (49) 



1 — V 2 dr prj dr 



4 An even better fit to the ideal relativistic Boltzmann gas equation of state is g(T) = (12 + 45T + 
45T 2 )/(8 + 2QT + 15T 2 ), which has a maximum relative error of 7 x 1CT 4 . 
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where 



1 Ail 9 - ■ 9 - 



f^-^'-nJ^-nJ (50) 



and f2 = vfi /u l = u + /P 1 ^ 2 /r 2 ^4 3//2 7, and Q± = ±(r 3 / 2 ± a) 1 is the rotation frequency, as 
observed at large radius, of circular, planar orbits. 

Our radial momentum equation obviously reduces correctly to the nonrelativistic limit. 
It agrees with PA's radial momentum equation, although theirs is expressed in different 
variables. It also agrees with ACGL's, although ACGL take rj — 1. 



3.4. Angular Momentum 

An angular momentum equation can be derived with the aid of the azimuthal Killing 
vector f£ = (0, 0, 0, 1). We have 

(T\$);» = = (T%), v = ±=.{JgT^), v = ^(r 2 T^), r + -^-{^n9T%) fi . (51) 

a / Q T lA olll U 

Vertically averaging, the second term vanishes, and integrating once and using the definition 
of the stress tensor, we obtain 

Mil] - A-nHerh^ = Mj. (52) 

Here Mj is the total inward flux of angular momentum; j is treated as an eigenvalue of the 
problem and must be self-consistently obtained. 

Our angular momentum equation agrees with PA and ACGL, except that ACGL 
assume rj — 1. 



3.5. Vertical Structure 

We assume vertical equilibrium, and use the expression for the simplified equilibrium 
scale height derived by ALP under the assumption that u 6 and u 6 t e are small. A convenient 
way to to express ALP's result is to define an effective vertical frequency via 



(53) 
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Then ALP find 5 

.iJl^A. ,54) 
In the limit of a relativistic, thin disk this reduces to v 2 z = r _3 (l — 4a/r 3 / 2 + 3a 2 /r 2 )/C, 



which agrees with the result of [Riffert fc Herold ( 1995 )| , and obviously reduces correctly to 
the nonrelativistic limit. 

ALP's result for the equilibrium scale height differs from PA's. The reason is that 
ALP neglect u 9 and its derivatives, while PA work in a quasi-cylindrical coordinate system 
(t, r, z, 4>) and thus neglect u z and its derivatives. The two assumptions are not equivalent. 
If the flow is quasi-spherical, then ALP's approach is superior. 

The major defect of all these approaches to vertical structure is that they assume 
vertical equilibrium. This assumption is not justified close to the event horizon. Inflow there 
is so rapid that inertial terms become important, and the scale height effectively "freezes 
out." ALP provide an expression for the evolution of the scale height that takes this effect 
into account, but we have not used this because it introduces significant mathematical 
complications. 



4. Turbulent Shear Stress 

The key to calculating the structure of steady, rotating accretion flows lies in 
understanding the nature of the viscous stress t r ' $ that allows the angular momentum of 
the flow to evolve. In our view the viscous stress is likely caused by MHD turbulence 
initiated by the Balbus-Hawley instability ( [Balbus fc Hawley ( 1991 )| ) . In this case, the 
"viscous" stress is really a time-averaged turbulent Reynolds and Maxwell stress (see 
Balbus, Gammic, fc Hawley (1994] ). 



The viscous stress is commonly assumed to have the Navier-Stokes form (see eq. 
One difficulty with this assumption is that, as pointed out by |Narayan (1992)| , the resulting 
system of equations are acausal: they propagate information at infinite speed. This has the 
undesirable consequence that one must apply a boundary condition at the event horizon 
(e.g. |Narayan, Kato, fc Honma (1997)|) . 



Various workers have sought to remedy this defect, including Narayan (1992), Narayan 



5 | This expression is closely related to a well-known integral of the motion for geodesies in the Kerr 
geometry: Q = u 2 g + cos 2 9 - a 2 (£ 2 - 1)) [Carter (1968)| . In the limit of small 69 = - tt/2, 



Q ~ Ug + r 4 88 2 v 2 . Evidently Q is proportional to a vertical energy of oscillation. 
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Loeb, & Kumar (1994), |Papaloizou fc Szuszkiewicz (1994")| (hereafter PS), [Kato fc Inagak: 



(1994)| , [Kato (1994)| , and Narayan (1996, private communication). We follow PS's line of 



reasoning because it is straightforward to generalize it to the relativistic case. 



4.1. Causal Viscosity: Nonrelativistic Case 

PS introduce a phenomenological equation for the evolution of the viscous stress that 
allows it to relax toward an equilibrium value. Let us review their result in our notation, 
temporarily assuming that the velocities are nonrelativistic and the geometry is Euclidean. 
Then PS's equation for the viscous stress 5* = f ' $ is 

DS S — S D . . 

Wr = — — ' (55) 

where D / Dr is the convective derivative, S Q is the equilibrium value of the stress and r r is 
the relaxation timescale. In a steady state 

v t = ~- < 56 > 

Now take d/dr of the angular momentum equation to obtain 

M dl ( „ dS\ , \ 



neglecting terms of order dHg/dr. Solve for dS/dr in equationflSTD, eliminate dS/dr from 
equation(|56|), and use S = — pvr 2 dVL / dr to find 

( r dn | V2 ldl ) m 
1 - 2rV/r I dr^clrdrh { } 



where c v = yv/r r is the speed at which viscous effects propagate (see the discussion of PS). 
If r r ~ (as one might expect for Balbus-Hawley induced turbulence), and v ~ ac^/fl, 
then c v ~ CgCx 1 / 2 . Rewriting equation(KSf) slightly, 



S = - pvr 



2rV/r 



idi f i __ Yl) 

r dr 1 c? y / r 2 



(59) 



Evidently this viscosity prescription introduces a singular point into the angular momentum 
equation at V 2 = c 2 ,, since the coefficient of dl/dr in the angular momentum equation 
vanishes there. This requires a new boundary condition at the singular point. The new 
boundary condition replaces the old, objectionable boundary condition at the surface of the 
accreting object. 
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4.2. Causal Viscosity: Relativistic Case 

We now generalize PS's causal stress prescription to the relativistic case. Begin 
with the viscous stress-energy tensor measured in the LRF, £( )(&)- We assume that all 
components vanish except t( r )(<£) — ^(<t>)(r)- This is the minimal modification of the inviscid 
equations that produces angular momentum transport. Notice that the LRF is the only 
safe place to make such modifications, since an arbitrary modification of the stress tensor in 
another frame can have unintended side effects. Then, with the aid of equations ( p6|) and 
(|8|), we find | 

= ^r\ r ADt ( , r) ^ ) = Ft {r)W . (60) 

We now identify £( r )(<£) with S, the component of the stress that obeys the relativistic 
analogue of equation (|55|): 

DS S — S D . . 

D^ = —^r ' (61) 

where now D/Dt = u^Q.^. Then in a steady state 

dS S — S . . 

u r — = . (62) 

dr r r 

The rest of the argument proceeds as before: we differentiate the angular momentum 
equation, solve for dS/dr, eliminate dS/dr from equation (1621), and solve for S. We find 



l-T r u r (2/r + d\nF/dr)\° F dr J ' 1 ' 

We must now specify S a . 

The most natural form for the equilibrium value of the LRF shear stress is the 
Navier-Stokes value: 

S = -2pr]ua {r){(j , h (64) 

where v = ac s rH 9 , following the usual a prescription. The calculation of ctmu) is a lengthy, 
but important, matter that is left to the Appendix. The result is also complicated and left 
in the Appendix. Here we abbreviate it by 

dl dV 

a (rm =a = a N + a L — + a v — . (65) 



We disagree with Abramowicz fc Prasanna (1990)| , who argue that the viscous stress (indeed, the 



viscosity) vanishes at the horizon. Measured in the LRF, the viscous stress will not generally vanish 
at the horizon because the horizon is not in any way a special location from the point of view of 
the fluid, which falls freely through it. Transforming to the BLF, one finds that for general t( a )(6)i 
V ' $ = 77V ' AD{t( r ){4>) + ftr@4>lrt(r){r))- F° r geodesies 7 ~ j r ~ I? -1 / 2 near the horizon, so if i( r )(<£) an d 
i(r)(r) are finite at the horizon, so is V ' 
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In the nonrelativistic limit, a = (l/2)rdQ/dr, while in the thin disk limit 
a = -(3/4)r- 3 / 2 £>/C ( [Novikov fc Thorne (1973)1) . Q 

Now we recover the full expression for the shear stress 5*. Substituting in equation (^3| 
for S , we find 



l-r r u r (2/r + dlnF/dr)\dr[ L+ eg fJ * F cgF dr 



(66) 



The coefficient of dl/dr vanishes, and so there is a singular point in the flow, where 



f^) 2 « cl 



4.3. Dissipation Function 

The dissipation function $ is best calculated in the LRF, where 

$ = -&^a {am (67) 

Since most of the components of t vanish in the LRF, and the metric is Minkowski there, 
one finds 



$ = -2Sa 



V c*vf dr J 

If a < and d(lrj)/dr > then $ is positive. 

It turns out that the dissipation function is not always positive in our solutions. In 
particular, a sometimes changes sign close to the horizon, an effect first noticed by |Anderson 



fc Lemos ( 1988 )|. Since the stress is not proportional to the shear rate, it need not change 



sign at the same point. The dissipation function then goes negative, indicating that the 
stress is transporting angular momentum outwards against the shear. In part this reflects 
the inadequacy of our phenomenological equation for the evolution of the viscous stress. In 
part this may be a physical effect, however, in the sense that small scale magnetic fields 
embedded in the flow are "wound up" by the shear and briefly unwind, reducing the internal 



7 One can show that — a is the maximum growth rate for the Balbus-Hawley instability in a thin Kerr 
disk. The factor D/C varies from 1 at large radius to 4/3 at the last stable orbit. In the mechanical analogy 
for the Balbus-Hawley instability developed by Balbus & Hawley (1992) this result is unaffected by radial 
motions, to first order in the eccentricity. 



-16- 



energy of the fluid (into which we folded the magnetic energy). But it seems unlikely that 
this unwinding is fully reversible. As we shall see in the sample numerical solutions, this 
effect has no practical consequences since it only occurs close to the horizon, where the 
inflow time is short compared to the heating time and so the flow is nearly adiabatic. 



5. Critical Points 

The accretion flow must pass through two critical points. The first is the usual sonic 
point, which is made rather complicated by the presence of relativistic terms. The second is 
the "viscous point" that originates in the vanishing of the coefficient of dl / dr in the angular 
momentum equation. Here we develop explicit expressions for the associated boundary 
conditions. 

Before proceeding, it is convenient to make two more approximations. These 
approximations eliminate terms that weakly couple the angular momentum equation to 
the other equations and enormously complicate the evaluation of the critical points. First, 
we replace dlaF/dr in equation fl6"6] ) by 1/ {r^\/V). This approximation is good to about 
15% inside r = 10 and far better outside. Second, we replace dlnrj/dr in equation (|66 ) 



by — (rj — l)/r]r. This approximation is good to about a factor of two. These terms 
are important mainly in the relativistic regime, where the shear stress is dynamically 
unimportant and the flow is nearly adiabatic. 



5.1. Sonic Point 

In order to obtain the sonic point conditions, we need to find an expression for the 
radial velocity gradient which involves no other derivatives. Therefore, we must find an 
explicit expression for the pressure derivative in ( fl9|) . To accomplish this, we first gather 
together the conservation equations for particle number and energy (|3^) and (|47|), and the 
vertical equilibrium equation (1531), and write them in differential form: 



dlnp + dlnHg + 2 + 1 dlnV + 1 - a 2 /r _ Q 
dr dr r 1 — V 2 dr r 2 T> 

W) -l)^l-^£ = J* (70) 
dr dr pu r 1 

2^ = (1 - Th(T)/rj)^ - 2 - - (71) 
dr dr r dr 
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We have used h(T) = 1 + g{T) + Tdg/dT, so that dr]/dr = h(T)dT/dr, and 
(1/ p)du/dT = h(T) — 1. Notice that we have retained a term proportional to dHo/dr; this 
term is assumed small, but we must retain it for consistency if we take equation Q56"D as our 
basic particle number conservation equation. 

Using (fTID to eliminate dlnHg/dr from (|36"D , we can solve for d In p/dr and <i In T/dr, 
and obtain dlnP/dr = dlnp/dr + dlnT/dr. Substituting the result into (fl9| ), we obtain 



y 2 dlnV 
1 - V 2 dr 

where 



T 



fi(T) 



' 1 d In z/ 2 
2 dr 



dlnV 



1 



1 - V 2 dr 



/ 3 (T) 



r 2 X> 



+ / 2 (T) 



1 - (V 2 )!! - Th(T)/i]) 



(72) 



• (73) 



/l(T) fe(T) - 1 + (1/2) (1 - Th(T)/ V ) ' h(T) - 1 + (1/2)(1 - T/i(T)/t7) 

In the nonrelativistic (T — > 0) limit, f\ = 270/(70 + 1) and fi = (70 — 1)/ (70 + 1). 

Finally, we must move all terms proportional to d In V/ dr from the right-hand to the 
left-hand side of this equation. This is more involved that it first appears because both 
$ and d In v 2 /dr include terms proportional to d In V/dr. First, $ can be simplified by 
substituting the shear stress S from fl5"2"D into fl68|), which allows us to write 

/$ 2af(lr] 



_ JJ_ 

pu r T FT 

We then divide a into a part which is proportional to dV/dr and a part which is not: 



a 



'df 
dr 



dV 
dr 



'df 
dr 



(74) 



(75) 



N 



where we have taken dl/dr = (dl/dr) N + (dl/dr)vdV/dr. The terms (dl/dr)^ and (dl/dr)v 
are derived from (B3) using the shear stress (|66|) . Thus, the din V/dr part of the $ term is 



h{T) 



cry + &L 



'dl\ \ 2f(lrj-j)VdlnV 
dr I „ / F dr 



Following a similar procedure, we find that the d In v 2 z jdr term has a dlnV/dr part 



(76) 



Tfi(T) (d\nvl 



2t] 



dr 



dV 
dr 



Tfi(T)V 



l-Sa 2 Q fdV 
r 4 \ dr 



Sa 2 (d£ 
r 4 \ dr 



v. 



d\nV 

dr 



(77) 



Solving for dlnV/dr, we obtain 



dlnV 

dr 



N. 



{7i 
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where 

N s = fr- 
aud 



Tfi (T) 



1 ( dhu/[ 

2 I dr 



iV 



r — 1 
r 2 X? 



/2(T) L + „ K 
77 V \dr 



(79) 



2VV 



^ = V 2 | T/x(T) fV_ (dlnvf 



1-V 2 



V 



dr 



1-V< 



f*(T) 



dP 

avJraL \Tr 



F 

(80) 

At the sonic point N s = D s = 0. In the nonrelativistic limit the condition D s = reduces 
to V 2 = c^27o/(7o + 1), in agreement with NKH. 



5.2. Viscous Point 

The viscous point conditions are obtained by solving the angular momentum equation 
(B2|) for dl/dr. We obtain 

dl _N V 

Tt-d: (81) 

where 

[l-r^/r + dlnF/dr)]^ 1 / 2 o ^ dV (u r ) 2 dln V 

Nv = F W (1 - VW 2 N V dR ~ -^F^T- (82) 

(u r ) 2 

D v = -2a L + (83) 
ctF 



and 



At the viscous point N v = D v = 0. The conditions !)„ = is equivalent to 



Yl 

c 2 



j(A/V) 1/2 (l-V 2 )(£-in). (84) 



We take c 2 = ac 2 , and since a < 1 the viscous point will generally be located well outside 
the sonic point. 



6. Example Numerical Solution 

In this section we find a typical numerical solution and use it to study the importance 
of various relativistic terms, to check our approximations, and evaluate the robustness of 
the input physics. More precisely, we have solved the system of equations fl59"|) , (|4"7|), (|4"5| ) 
and (|5^) numerically, using a relaxation scheme (see [Popham fc Gammie (1997)1 f° r more 
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details). We have set / = 1 (zero cooling), a = 0.1, and the rotation parameter of the black 
hole a = 0. 

The system is subject to five boundary conditions that allow us to solve the first order 
differential equations for V, I, p, T and for the eigenvalue j over a domain that extends 
from just outside the event horizon to r = 2 x 10 A GM/c 2 . Two boundary conditions are 
provided by equations (|79|) and (j82|) , which ensure the solution passes smoothly through 
the sonic and viscous points, respectively. Two additional boundary conditions fix I and T 
at the outer edge of the solution (V must be determined self-consistently there). We set 
T ut = 1-67 x 10~ 5 and l ou t = 57.75. These values are taken from the self-similar solution 
of |Narayan fc Yi (1994)| . NKH have shown that global solutions with, e.g., thin disk outer 
boundary conditions, rapidly approach this self-similar solution, which then spans many 
decades in radius. Consistent with this, we find that the character of the solution at r ^ 10 2 
is insensitive to the outer boundary conditions. A final boundary condition is obtained by 
normalizing the density so that the rest-mass accretion rate M — 1. 

The run of the basic variables V, I, p, T with radius is shown in Figure 1. The upper 
left panel shows V, and the arrows indicate the location of the sonic point at r s = 6.41 
(slightly outside the last stable circular orbit at r = 6) and the viscous point at r v = 28.2. 
At the viscous point V ~ 0.3c s and so it lies far outside the sonic point. The solution is, 
however, insensitive to the precise definition of c„. For example, setting c u = 0.5c s changes 
j by less than 4%. The solution reaches V = 1 at the event horizon, which is located at 
r = 2. 

The upper right panel in Figure 1 shows I (heavy solid line), Irj (dashed line), and the 
eigenvalue j = 2.62 (dotted line). Recall that It], rather than /, is conserved in the absence of 
viscous torques, and that Mj is the total inward flux of angular momentum. The difference 
If] — j is proportional to the shear stress; in this case the difference is a small positive 
number at the horizon. The lower left panel shows density, which very nearly follows a 
power law r~ 3 / 2 over the whole solution. The lower right panel shows the temperature, 
which reaches a maximum of 0.067, corresponding to a sound speed c s = 0.28c. As we shall 
see, the maximum temperature and density of the solution change significantly with black 
hole rotation. 

The run of some important dynamical quantities is shown in Figure 2. The upper left 
panel shows the absolute value of the shear rate in the LRF. The shear changes sign at 
r w 3. Also shown is the thin disk approximation to the shear rate, equation (|A6|) , which 



8 Physical units may be recovered as follows: radial velocity is Vc, angular momentum is IGM/c, density 
is pMG/c 3 (since the mass accretion rate is normalized to 1), and temperature is Tine 2 /k. 
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lies too close to the actual shear rate to distinguish it on this plot. The maximum difference 
between the two is about one part in 10 3 . The upper right panel shows the causal shear 
stress as a solid line and the equilibrium shear stress S (proportional to the shear rate, 
see equation p3J) as a dashed line. The causal stress is relaxing toward So but the fluid is 
also flowing inwards toward regions of larger shear, so the causal stress always lies slightly 
below S*o. The lower left panel shows Hq. The vertical averaging approximation is accurate 
to of order dlnHg/dlnR, which is largest at the horizon. The lower right panel shows 
Q = u^/u 1 . For a = we must have fl = at the horizon, so Q goes through a maximum at 
r m 3. This is consistent with our statement that the shear rate in the LRF is very nearly 
proportional to dQ/dr. 

Thermodynamic quantities are shown in Figure 3. These are the relativistic enthalpy 
rj, the pressure p, the entropy S (minus its value at the outer edge), and the dissipation 
function $. As discussed earlier, the dissipation function goes negative at r w 3.0, where 
the shear changes sign. This is an inevitable result of our causal viscosity prescription. It 
has little practical significance, however, since, as can be seen from Figure 3, the flow is 
nearly isentropic at r < 3.0. This is because the radial inflow time becomes short compared 
to the heating time. 

This work contains two new physical ingredients as compared to NKH: a causal 
viscosity prescription, and a full treatment of relativistic effects. We can evaluate the 
influence of each by solving a simplified version of our equations which retains our causal 
viscosity formulation, but eliminates all relativistic corrections and terms and adopts the 
Paczyhski potential used by NKH. We can compare the NKH solution and the causal 
solution to find the effects of the causal viscosity, and we can compare the causal solution 
and the full, relativistic solution to find the effect of properly including relativity. 

The results are displayed in Figure 4, which shows the basic variables V,l,p,T in 
each of the three solutions: the fully relativistic solution (heavy solid line), the causal 
solution (long dashed line) and the NKH solution (dotted line). Also shown in the angular 
momentum plot is a short dashed line, which is the angular momentum of a thin disk 
(inside r = 6 the angular momentum is that of the last stable orbit). Before continuing 
with the comparison, notice that all the ADAF solutions lie below the thin disk solution in 
angular momentum- evidently they are all "sub-Keplerian" . 

Remarkably, the relativistic solution, causal solution, and NKH solution are very 
similar at r ^ 10. This suggests that the solutions are robust in the sense that significant 
changes in the input physics cause only small changes in the solution. One can also check 
the robustness of the shear stress model by changing the definition of c 2 ,. If we set c 2 = c 2 
(rather than ac 2 s ) we find that the values of the basic dependent variables change by no 
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more than 5.9% at the inner boundary. Again, this suggests that the solution is robust. 

What is the effect of the causal viscosity? The comparison of the causal solution 
and the NKH equations is not perfect because NKH take v\ = Q 2 K (which diverges at 
r = 2), while we take v\ = Q 2 . The outcome of combining these two changes is to lower 
the angular momentum slightly and decrease the radial velocity (therefore increasing the 
density) in comparison to NKH. What is the effect of properly including the relativity? The 
most obvious effect is in the radial velocity, which diverges at r = 2 in the causal solution 
because the potential gradient diverges there. This causes the density and temperature to 
go through a maximum and then decline to zero at r = 2. In the relativistic solution, by 
contrast, the density and temperature increase monotonically in to the horizon. The next 
most obvious effect is in the angular momentum profile, which is significantly lower at the 
horizon in the relativistic solution than in the causal solution. This is due to our retention 
of relativistic enthalpy terms in the angular momentum equation; the angular momentum 
per unit rest mass is Irj rather than I. 

The similarity of all three ADAF solutions in Figure 4 should not lead one to conclude 
that relativistic effects are not important. Figure 5 shows the effect of varying a from 
—0.999 to to 0.999. The temperature at the horizon increases from 0.067 when a = to 
0.363 when a = 0.999, while the density increases from 0.11 at a = to 1.29 at a = 0.999. 
These changes suggest the possibility that ADAF emission properties may depend strongly 
on black hole rotation. We defer a more detailed discussion of the effect of black hole 
rotation to Popham fc Gammie (1997) . 



We can now check our approximations for self-consistency. Vertical averaging cannot 
be checked directly, but the smallness of | d In Hg/d In r| compared to 1 is a rough test of 
its validity (this also tests the validity of our vertical equilibrium approximation). The 
maximum in \d In Hg/d In r| is 0.697 at the horizon and much smaller outside (see Figure 2). 
The situation would be improved if the scale height were allowed to evolve dynamically, so 
that Hq could "freeze out" near the horizon. 

We have also dropped viscous accelerations in the radial momentum equation. One 
would expect that the viscous term is of order a compared to the dominant terms, and 
a direct check confirms this. The viscous term generally produces an outward directed 
acceleration, although in the inner parts of those solutions where the shear changes sign 
(generally when a ^ 0.7) it produces an inward acceleration. 

We made two additional approximations prior to deriving the sonic point conditions: 
(1) we set dlnF/dr ~ l/(r7\/P); (2) d\nr]/dr ~ — (rj — l)/r]r. We have checked these 
approximations by multiplying the approximated term by 0.5 and examining the effect on 
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the solution. For approximation (1) we find the largest fractional change is 2.7% in T at the 
innermost radius. For (2) we find the largest fractional change is 0.2% in / at the innermost 
radius. This shows that the approximations, whose accuracy was discussed at the beginning 
of §5, do not affect the solution very much. 

Finally, how important are the relativistic thermodynamic terms (rj and g) that we 
have worked so hard to include? The rj term turns out to be quite important, particularly 
for holes with rapid prograde rotation. When a = 0.999, the maximum value of rj is 
2.26; this implies substantial corrections in the angular momentum and radial momentum 
equations. The relativistic equation of state is somewhat less important; the maximum 
ratio of g to 1/ (7 C — 1) (the nonrelativistic value) is 1.10. 



7. Summary 

We have written down and solved equations describing optically thin, advection 
dominated flows in the Kerr metric. Our physical description of the flow includes a causal 
prescription for the viscosity, and does not require the application of a boundary condition 
at the event horizon. We solve the energy equation assuming that a constant fraction of the 
dissipated energy is advected. We have also included certain relativistic terms that have 
been neglected in some earlier treatments; most importantly, we do not assume that the 
relativistic enthalpy rj is unity. We will present a detailed survey of solutions for various 
choices of parameters in an upcoming paper ( |Popham fc Gammie (1997)| ). These solutions, 



when coupled with a detailed description of the cooling processes and a photon transport 
scheme, will allow us to produce model spectra for advection dominated systems. 



We are grateful to R. Narayan for helping to initiate this project and for his support, 
insight and encouragement, and to J.-P. Lasota and W. Press for helpful discussions. This 
work was supported by NASA grant NAG 5-2837 and NSF grant AST 9423209. 



A. Shear Rate in the Local Rest Frame 



We want vmu), which is the shear rate measured in the local rest frame (LRF). Using 
the basis vectors for the LRF, 

°"(r)(i£) = e^e^cr^, (Al) 
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where is given by equation (|33|). We must then calculate many of the pieces of a^, 
although our task is simplified a bit because 



e /1 (a)^ = 5{*)> 



(A2) 



since e^u) is parallel to it^, and the basis vectors are orthonormal. Along the way we find 



u t -t 

Ut-r 
U t -4, 



dS 

dr 

a(3 r ~j r \/V 



7 A 



du r 
dr 



+ /3 r7r (l -a 2 /r) 



-1 

Ut-4, 
dl 



a-a 2 Q + r 3 Q) 



rf3 r7r W(l-a 2 /r 3 ). 



Using equation (JT 

dS 

dr 



VdV dl luo 



+ n—-^-(3 + a 2 /r 2 ) + 



7(1 — a 2 /r) 



7 V A dr dr rA r 
7v / Pa 2 (l + 3/r) / 2 (1 - a 2 /r 3 ) F5 

M 



+ 



3^43/2 



r 3 A 2/ y 



Assembling the results and transforming to the LRF, we find 



IV 1± (VdV E-lSldl l 2 oo(3 + a 2 /r 2 ) 1 2 , , . 2N 

7^(3 + a 2 /r 2 ) 



+ 



+ 



2r7 V -4. cir 2r dr 2r 2 ^ 2 A 2 

l£[An_V(\- a 2 /r 3 ) 
7 r V ^ \ r3 



cu 2 (2 + a 2 /r 2 + a 2 /r 3 ) 



(A3) 



(A4) 



(A5) 



| l 2 oo{3 + a 2 /r 2 ) / 7 (l-a 2 /r) / 7 y / Pa 2 (l + 3/r) | Z 3 (l - a 2 /r 3 ) /P 



2r 2 ^l 2r 3 v^4P 2rM 3 / 2 2rM 2 7 V A' 

For a relativistic thin disk, this reduces to the remarkably simple expression 



1 . 2 dVL 



(THIN DISK) 



(A6) 
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( |Novikov fc Thorne (1973] ). Obviously this reduces correctly to the nonrelativistic limit. 



Remarkably, this thin disk expression for a lies within a fraction of a percent of the full 
LRF shear stress in our advection dominated solutions. 

For cold, geodesic flow (/ = const., £ = const.) onto a nonrotating hole, a reduces to 



a (rm = -^=-r(l - 3/r). (a = 0, GEODESIC) (A7) 

This case is simple enough to be checked by hand. It shows the shear reversal at r = 3 
discovered by |Anderson fc Lemos (1988)1 an d discussed in detail by |Abramowicz fc Prasanna 

rrggm 
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10° 10 1 10 2 10 3 10 4 10° 10 1 10 2 10 3 10 4 

r c 2 /G M r c 2 /G M 

Fig. 1. — Solution to the full relativistic equations with a = 0, a = 0.1, / = 1. The upper 
left panel shows the radial velocity V measured in the CRF. The arrows show the location of 
the sonic and viscous critical points. The upper right panel shows the angular momentum I 
(solid line) as well as Irj (dashed line) and j (dotted line). The lower left panel shows density 
and the lower right panel shows dimensionless temperature T. 
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Fig. 2. — Dynamical quantities in the a = solution. The upper left panel shows the shear 
rate measured in the LRF. The approximate formula for the shear rate, based on the thin 
disk shear, is also shown but is too close to distinguish on this plot. The upper right panel 
shows both the causal shear stress (solid line) and the acausal (Navier-Stokes) version, which 
is proportional to the shear rate (dashed line). The lower left panel shows the angular scale 
of the flow He, while the lower right shows the angular frequency Q = /u*. 



-29- 




10° 10 1 10 2 10 3 10 4 
r c 2 /G M 



lO" 2 


J_ 1 1 | Mill 


Ilj 1 1 llllllj 1 1 llllll 


4 


io- 3 








io- 4 






~i 


io- 5 


r \ 






io- 6 








io- 7 


r 




-1 


io- 8 








io- 9 








10 -io 


r 






io- 11 










S 1 


III 1 


1 if 



10° 10 1 IO 2 10 3 10 4 



r c 2 /G M 



io- 3 


jj- 1 1 1 lllllj — 1 1 1 lllllj — 


1 1 llllllj 1 1 1 lllllj TF| 


10- 4 




n 


io- 5 


r |\ 


"i 


io- 6 






IO" 7 




-| 


IO" 8 






io- 9 






*r 10- 10 
— 10" 11 


r 




io- 12 


r 


\ H 


io- 13 


r 


\ n 


10 -14 






IO" 15 






10 -16 


r 




io- 17 


» I I 


1 1 i i f 



10° 10 1 10 2 10 3 10 4 

r c 2 /G M 



Fig. 3. — Thermodynamic quantities in the a = solution. The upper left panel shows the 
relativistic enthalpy rj, the upper right panel shows pressure, the lower left entropy, and the 
lower right the absolute value of the dissipation function. The dissipation function changes 
sign at r ~ 3. 
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Figure 4 
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Fig. 4. — A comparison of various ADAF solutions: the relativistic solution (heavy solid 
line); a solution of the causal nonrelativistic limit of the relativistic equations (long dashed 
line); the NKH solution (dotted line); and the thin disk angular momentum distribution 
(short dashed line). 
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Fig. 5. — A comparison of the relativistic solution for different black hole angular momenta. 
Shown are solutions for a = (solid line), a = —0.999 (long dashed line) and a = 0.999 
(short dashed line). 



